# substantive effects background
load("estimates.RData")
source("effectSim.R");source("effectSimPlot.R")
library("AER");library(arm)
summary(modslib2[[1]])
x <- cbind(1,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1:40,1:40,1:40,1:40,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
baseline <- effectSim(model=modslib2[[1]],covvar=vcovHAC(modslib2[[1]]),
                      x=x,sims=1000)
x1 <- cbind(1,0,2,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1:40,1:40,1:40,1:40,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0)
westfarm <- effectSim(model=modslib2[[1]],covvar=vcovHAC(modslib2[[1]]),
                      x=x1,sims=1000)
input.prob <- function(p,low,high){
  predicted.prob <- t(apply(p, 1, quantile, probs = c(0.05, 
                                                      0.5, 0.95)))
  cord.y <- c(predicted.prob[, 1], rev(predicted.prob[, 
                                                      3]))
  cord.x <- c(seq(low, high, length = nrow(p)), seq(high, 
                                                    low, length = nrow(p)))
  return(cbind(cord.x,cord.y))
}
west <- input.prob(westfarm,low=1,high=40) 

summary(modscon2[[1]])
x2 <- cbind(1,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1:40,1:40,1:40,1:40,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
baselineC <- effectSim(model=modscon2[[1]],covvar=vcovHAC(modscon2[[1]]),
                      x=x2,sims=1000)
x3 <- cbind(1,0,2,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,1:40,1:40,1:40,1:40,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
french.prob <- effectSim(model=modscon2[[1]],covvar=vcovHAC(modscon2[[1]]),
                      x=x3,sims=1000)
french <- input.prob(french.prob,low=1,high=40) 

lab <- c(1867,1872,1874,1878,1882,1887,1891,1896,1900,1904,1908,1911,1918,1922,1925,1926,1930,1935,1940,1945,1949,1953,1957,1958,1962,1963,1965,1968,1972,1974,1979,1980,1984,1988,1993,1997,2000,2004,2006,2008)
tlab <- 1:40

pdf("Figure-4a.pdf",width=7,height=2.75)
par(mar=c(4.1,4.1,3.1,2.1))
effectSimPlot(baseline,low=1,high=40,i=21,x=x,ylim=c(.3,1),rug=F,color="grey70",
              covariate=allout$parliament,xlab="Liberals"
              ,ylab="Voting Loyalty",bty="n",xaxt="n")
polygon(west[,1],west[,2],col="grey40",border=NA)
lines(x[, 21], t(apply(westfarm, 1, quantile, probs = 0.5)), lty = 1)
axis(1,at=seq(1,40,1),labels=lab)

legend("right",legend=c("Ontario","Western"),
       col=c("grey70","grey40"),lwd=4,bty="n",lty=c(1,1,3),cex=.75)
dev.off()
pdf("Figure-4b.pdf",width=7,height=2.75)
par(mar=c(4.1,4.1,3.1,2.1))
effectSimPlot(baselineC,low=1,high=40,i=21,x=x,ylim=c(.3,1),rug=F,color="grey70",
              covariate=allout$parliament,xlab="Conservatives"
              ,ylab="Voting Loyalty",bty="n",xaxt="n")
polygon(french[,1],french[,2],col="grey40",border=NA)
lines(x[, 21], t(apply(french.prob, 1, quantile, probs = 0.5)), lty = 1)

axis(1,at=seq(1,40,1),labels=lab)

legend("right",legend=c("English speaking (Ontario)","French speaking (Quebec)"),
       col=c("grey70","grey40"),lwd=4,bty="n",lty=c(1,1,3),cex=.75)
dev.off()
